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ABSTRACT 

Context. According to the sequential accretion model (or core nucleated accretion model), giant planet formation is based first on the 
formation of a solid core which, when massive enough, can gravitationally bind gas from the nebula to form the envelope. The most 
critical part of the model is the formation time of the core: in order to trigger the accretion of gas the core has to grow up to several 
Earth masses before the gas component of the protoplanetary disc dissipates. 

Aims. In this paper, we calculate planetary formation models including a detailed description of the dynamics of the planetesimal 
disc, taking into account both gas drag and excitation of forming planets. 

Methods. We compute the formation of planets, considering the oligarchic regime for the growth of the solid core. Embryos growing 
in the disc stir their neighbour planetesimals, exciting their relative velocities, which makes accretion more difficult. Here we introduce 
a more realistic treatment for the evolution of planetesimals' relative velocities, which directly impact on the formation timescale. For 
this, we compute the excitation state of planetesimals, as a result of stirring by forming planets, and gas-solid interactions 
Results. We find that the formation of giant planets is favoured by the accretion of small planetesimals, as their random velocities are 
more easily damped by the gas drag of the nebula. Moreover, the capture radius of a protoplanet with a (tiny) envelope is also larger for 
small planetesimals. However, planets migrate as a result of disc-planet angular momentum exchange, with important consequences 
for their survival: due to the slow growth of a protoplanet in the oligarchic regime, rapid inward type I migration has important 
implications on intermediate mass planets that have not started yet their runaway accretion phase of gas. Most of these planets are lost 
in the central star. Surviving planets have either masses below 10 M e or above several Jupiter masses. 

Conclusions. To form giant planets before the dissipation of the disc, small planetesimals (~ 0. 1 km) have to be the major contributors 
of the solid accretion process. However, the combination of oligarchic growth and fast inward migration leads to the absence of 
intermediate mass planets. Other processes must therefore be at work in order to explain the population of extrasolar planets presently 
known. 

Key words. Planets and satellites: formation - Planet-disc interactions - Methods: numerical 



1. Introduction 

Since the discovery of the first extrasolar planet around a solar- 
type star (Mayor & Queloz, 1995) about 800 extrasolar planets 
have been identified. Observations indicate that planets are abun- 
dant in the universe. Planets orbiting stars show a great variety 
of semi-major axis (from less than 0.01 AU to more than hun- 
dreds of AU) and masses (from less than an Earth mass to sev- 
eral Jupiter masses), as can be found in The Extrasolar Planets 
Encyclopaedia (http://exoplanet.eu/). Planet formation models 
should be able to explain this observed diversity. The sequential 
accretion model or also called core nucleated accretion model 
is currently the most accepted scenario for planetary formation 
(e.g. Mizuno 1980, Pollack et al. 1996, Alibert et al. 2005 a, b, 
among others), as it can account naturally for the formation of 
planets in all mass ranges^] It proposes that planetary growth 
occurs mainly in two stages. In the first stage, the formation 
of planets is dominated by the accretion of solids. If the pro- 
toplanet is able to grow massive enough (~ 10 M ffi ) while the 
gas component of the protoplanetary disc is still present, it can 
bind gravitationally some of the surrounding gas, giving birth to 
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1 Note, however, that the formation of planets at large distance from 
their central star seems very difficult in the core-accretion model. 



a gas giant planet. The accretion of gas is slow at the beginning: 
the planet growth is dominated by the accretion of solids and 
the energy released from the accreted planetesimals slows down 
the accretion of the envelope. When the accretion of planetesi- 
mals declines (generally because the protoplanet has emptied its 
feeding zone) the accretion of gas is triggered, and the planet 
can accrete hundreds of Earth masses in a very short time. The 
runaway accretion of gas characterises the second stage of the 
sequential accretion model, where the growth of the planet is 
dominated by the accretion of gas. 

Since it was first proposed (Mizuno 1980), the sequential ac- 
cretion model has been extensively studied and improved, trying 
to include the many fundamental processes that occur simulta- 
neously with the growth of the planet, and that impact directly 
on it. Constructing a complete model that accounts for all these 
processes in a reasonable way is a hard task. Among its main in- 
gredients, it has to include a realistic model for an evolving pro- 
toplanetary disc, a model for the accretion of solids and gas to 
form the planets (which itself requires a knowledge of the inter- 
nal structure of the planet), a model for the interactions between 
the planets and the disc, and a model for the interactions between 
the forming planets. Each of these topics represent itself an in- 
dependent, ongoing area of research. Alibert et al. (2005) (A05 
from now on) include some of these processes in a single planet 
formation model. Given the complexity of the problem or the 
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unknowns related to some of these processes, many simplifica- 
tions have to be assumed in order to keep the problem tractable 
from the physical and computational point of view. This is spe- 
cially important, as we aim to compute thousands of simulations 
in order to account for the wide range of possible initial condi- 
tions (see Mordasini et al. 2009, M09 from now on, for more 
details). Therefore, our models represent a compromise between 
accuracy and simplicity in their physical description. 

The first stage of planetary formation corresponds to the 
growth of the solid embryo, which is dominated by the accretion 
of planetesimals. The growth of an embryo proceeds itself in two 
different regimes (Ida & Makino 1993, Ormel et al. 2010). At 
the beginning, big planetesimals, that have larger cross-sections, 
are favoured to grow even bigger by accreting planetesimals that 
they encounter on their way. Being more massive, in turn, en- 
larges the gravitational focusing, which leads to accretion in a 
runaway fashion. However, at some point, these runaway em- 
bryos become massive enough to stir the planetesimals around 
them. This results in an increase of the relative velocities and the 
corresponding reduction of the gravitational focusing. Growth 
among small planetesimals is stalled and only big embryos have 
the possibility to continue accreting, although at a slower pace. 
This second regime in the growth of solid embryos is known as 
oligarchic growth, as only the larger planetesimals or embryos 
(the oligarchs) are able to keep on growing. One important as- 
pect is that the transition between runaway and oligarchic growth 
occurs for very small embryos. As shown in Ormel et al. (2010), 
the actual mass for this transition depends upon many factors: 
size of the accreted planetesimals, location in the disc, surface 
density of solids, among others. In most of the cases, an em- 
bryo of ~ 0.01 M ffi , or even smaller, is already growing in the 
oligarchic regime. 

Our model is building up on the model of A05 and M09. 
Our primary aim is to study the formation of planets of differ- 
ent sizes, and in particular the cases for which the accretion of 
gas is important. Therefore, in the computations presented here, 
we focus on the first phase of planetary formation, when the gas 
component of the disc is still present (in the case of small rocky 
planets, collisions between embryos after the dissipation of the 
disc should be included in order to calculate their final masses). 
For the formation of giant planets, the growth of the solid core 
is dominated by the oligarchic growth. One of the weak points 
of the majority of previous giant planet formation models (e.g. 
Pollack et al. 1996, Hubickyj et al. 2005, A05, M09, Lissauer 
2009, Mordasini et al. 2012, hereafter M12) is the description 
of the solid disc, in particular of the interactions between form- 
ing planets and planetesimals. Such simplified models lead to an 
overestimation of the solid accretion rate which, in turn, results 
in an underestimation of the formation time of the whole planet. 

Indeed, in those works, the model for the accretion rate 
of solids is oversimplified: the whole formation of giant plan- 
ets' cores is assumed to proceed very fast, underestimating the 
excitation that planetesimals suffer due to the presence of the 
embryos. When oligarchic growth is adopted as the dominant 
growth model, giant planet formation turns out to be more diffi- 
cult. Formation times become much longer than the typical life- 
time of the protoplanetary disc. Fortier et al. (2007, 2009) and 
Benvenuto et al. (2009) studied the formation of giant planets 
adopting the oligarchic growth for the core. Assuming in situ 
formation for the planets and a simple, non evolving protoplan- 
etary disc, they show that the formation of giant planets is un- 
likely if planetesimals populating the disc are big (more than a 
few kilometres in size). However, formation could be acceler- 
ated if most of the accreted mass is in small planetesimals (less 



than 0.1 km). Guilera et al. (2010, 2011), also considering in situ 
models, studied the simultaneous formation of several planets 
where planetesimal drifting is included. They consider different 
density profiles for the disc and they find that, only in the cases 
of massive discs, the formation of the giant planets of the Solar 
System is possible only if planetesimal radii are smaller than 1 
km. These models, however, do not take into account that planets 
would likely migrate during their formation. 

In this work we include in our planet formation model a more 
realistic description for the accretion of solids. In Sect. [2] we re- 
view the basics of the A05 formation model, presenting some 
improvements in the computation of the disc structure, internal 
structure and migration. In Sect. [3] we describe the new treat- 
ment of the accretion of planetesimals. In Sect. |4] we present the 
results obtained for the formation of isolated planets (the forma- 
tion of planetary systems is described in Alibert et al. 2012, and 
Carron et al. 2012). In Sect. [5] we discuss our results and put 
them in context. Finally, in Sect.[6]we summarise our results and 
underline the main conclusions. 



2. Formation model 

The model and the numerical code used to calculate the forma- 
tion of planets is in essence the same as in A05. In what fol- 
lows, we make a summary of the most relevant aspects of the 
model and the improvements that have been introduced since 
that work. In the next section, we focus on the accretion rate 
of solids and we describe in detail the adopted model for the 
protoplanet-planetesimals interactions. 



2. 1. Protoplanetary disc: gas phase 

The structure and evolution of the protoplanetary disc is com- 
puted by first determining the vertical structure of the disc, for 
each distance to the central star, and second, computing the ra- 
dial evolution due to viscosity, photoevaporation, and mass ac- 
cretion by forming planets. 



2.1.1. Vertical structure 

The vertical disc structure is computed by solving the following 
equations: 



1 8P 

Pgas dl 

dF 9 
~dz ~ 4 



= -&z, 



p gas vfi 



and 
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3k P . 



gas 



dz 
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They reflect the hydrostatic equilibrium, the energy conserva- 
tion, and the diffusion for the radiative flux. In these equations, 
z is the vertical coordinate, p gas the gas density, P the pressure, 
T the temperature, v the macroscopic viscosity, F the radiative 
flux, k is the opacity (Bell & Lin 1994), and cr SB is the Stefan- 
Boltzmann constant. The Keplerian frequency, £2, is given by 



Q 2 = GM*/a 3 , 



(4) 
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G being the gravitational constant, M* the mass of the central 
star and a the distance to the stat0 

The equations ([T])-(|3]l are solved with four boundary condi- 
tions. The first three are the temperature, the pressure and the en- 
ergy flux at the surface. The surface of the disc is defined as the 
place where the vertical optical depth (between the surface and 
infinity) is equal to 0.01. The fourth boundary condition is that 
the energy flux equals in the midplane (see A05 for details). 
The three differential equations, together with the four boundary 
conditions, have a solution only for one value of the disc thick- 
ness H which gives the location of the disc surface. The macro- 
scopic viscosity v is calculated using the standard Shakura & 
Sunyaev (1973) a-parametrization, v = acJ/Q. The speed of 
sound c s is determined from the equation of state (Saumon et 
al. 1995). The temperature at the surface r sur f is computed as in 
A05. In the models presented in this paper, a is set to 7xl0~ 3 . 
This value of the alpha parameter has to be taken as an example. 

In the calculations of this work we neglect irradiation and the 
possible presence of a dead zone. These effects will be included 
in future works. 



2.1.2. Evolution 

rH 

The evolution of the gas disc surface density (S = J_ w p gas afz) is 
computed by solving the diffusion equation: 



dt 



3d_ 
a da 



da 



+ t w {a) + QplanetO). 



(5) 



where v is the effective viscosity, v = ^ ^ vp gas t/z . 
Photoevaporation is included using the model of Veras & 
Armitage (2004): 



t, w = for a < R g , 
S w oc or 1 for a > R„, 



(6) 



where R g — 5 AU. The total mass loss due to photoevaporation 
is a free parameter. The sink term Qpianet 1S equal to the gas mass 
accreted by the forming planets. For every forming planet, mass 
is removed from the protoplanetary disc in an annulus centred 
on the planet, with a width equal to the planet's Hill radius 



- a M 



\3Mj 



1/3 



(7) 



where M is the total mass of the planet and om is the location of 
the planet. 

Eq. [5] is solved on a grid which extends from the innermost 
radius of the disc to 1000 AU. At these two points, the surface 
density is constantly equal to 0. The innermost radius of the disc 
is of the order of 0.1 AU. 

Fig.[T]presents a typical evolution of a disc, whose parame- 
ters correspond to the first row of table [T] where the curves are 
plotted every 10 5 years. In this model, the photoevaporation term 
is adjusted in order to obtain a disc lifetime equal to 3 Myr. 

The characteristics of the protoplanetary disc are chosen to 
match as close as possible the observations. The initial disc den- 
sity profiles we consider are given by: 



£ = (2 - y) 



M di; 
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2 Note that we assume that the disc is thin, and the distance to the 
central star does not vary on a vertical slide of the disc. 




Fig. 1. Disc model. Gas surface density as a function of radius at 
different times. For the initial model (red line) the parameters are 
those of the first row of Table [T] the disc lifetime being 3 Myr. 
Each line corresponds to a different time (from top to bottom 
to 3 Myr, every 0.1 Myr). 

Table 1. Characteristics of disc models 



disc 


M disc (M Q ) 


a c (AU) 


dinner (AU) 


7 


1 


0.029 


46 


0.14 


0.9 


2 


0.117 


127 


0.16 


0.9 


3 


0.143 


198 


0.10 


0.7 


4 


0.028 


126 


0.10 


0.4 


5 


0.136 


80 


0.10 


0.9 


6 


0.077 


153 


0.12 


1.0 


7 


0.029 


33 


0.10 


0.8 


8 


0.004 


20 


0.10 


0.8 


9 


0.012 


26 


0.10 


1.0 


10 


0.007 


26 


0.10 


1.1 


11 


0.007 


38 


0.10 


1.1 


12 


0.011 


14 


0.10 


0.8 



where a ( ) is equal to 5.2 AU. The mass of the disc (M^ sc ), the 
characteristic scaling radius (ac) and the power index (y) are de- 
rived from the observations of Andrews et al. (2010). Adopting 
this kind of initial density profile is a difference from previ- 
ous works (A05 and M09). For numerical reasons, the inner- 
most disc radius, fli nne r, is always greater than or equal to 0.1 
AU, and differs in some cases from the one cited in Andrews 
et al. (2010). The afore-mentioned parameters used to generate 
the initial disc's profile are listed in table [T] Note that Andrews 
et al. (2010) derive also a value for the viscosity parameter or. 
On the contrary, and for simplicity, we assume here that the vis- 
cosity parameter is the same for all the protoplanetary discs we 
consider (a — 7 x 10~ 3 ). Using a different a parameter will be 
the subject of future work. Note also that, in the observations of 
Andrews et al. (2010), the mass of the central star ranges from 
0.3 to 1.3 M Q . However, we assume here that these disc pro- 
files are all suitable for protoplanetary discs around solar mass 
stars. Future disc observations will help improving this part of 
our models. 

As in A05, the planetesimal-to-gas ratio is assumed to scale 
with the metallicity of the central star. For every protoplanetary 
disc we consider, we select at random the metallicity of a star 
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from a list of ~ 1000 CORALIE targets (Santos, private com- 
munication). 

Finally, following Mamajek (2009), we assume that the cu- 
mulative distribution of disc lifetimes decays exponentially with 
a characteristic time of 2.5 Myr. When a lifetime T^ sc is selected, 
we adjust the photoevaporation rate in order that the protoplan- 
etary disc mass reaches 10~ 5 Mq at the time t = T^ c , when we 
stop the calculation. 

2.2. Protoplanetary disc: solid phase 

We consider that the planetesimal disc is composed of rocky and 
icy planetesimals. Here we assume a mean density of 3.2 g/cm 3 
for the rocky ones, and 1 g/cm 3 for the icy. The rocky planetes- 
imals are located between the innermost point of the disc (given 
by the fourth column of table [TJ, and the initial location of the 
ice line, whereas the disc of icy planetesimals extends from the 
ice line to the outermost point in the simulation disc. 

The location of the ice line is computed from the initial gas 
disc model, using the central temperature and pressure. The ice 
sublimation temperature we use depends upon the pressure. Note 
that in our model, the location of the ice line does not evolve with 
time. In particular, no condensation of moist gas, or sublimation 
of icy planetesimals is taken into account. Moreover, the location 
of the ice line is based on the central pressure and temperature, 
meaning that the ice line is taken to be independent of the height 
in the disc. Note that in reality the ice line is likely to be an "ice 
surface" whose location depends upon the height inside the disc 
(see Min et al. 2011). 

For the models presented here, we assume that all planetesi- 
mals have the same radius. Planetesimals' mass is calculated as- 
suming that they are spherical and have constant density (which 
depends on their location in the disc), and it does not evolve with 
time. The extension of our calculations towards non-uniform and 
time evolving planetesimal mass function is something we are 
working on and will be included in the next paper. 

The surface density of planetesimals, E m , is assumed to be 
proportional to the initial surface density of gas This means 
that 



^m(fl) = /D/G/R/i(a)£o(a), 



(9) 



where /d/g is the dust-to-gas ratio of the disc and it scales with 
the metallicity of the central star (the lowest value being 0.003 
and the largest one 0.125), and /r/i takes into account the de- 
gree of condensation of heavy elements. As in M09, we consider 
/r/i = 1/4 inside the ice line, and f R /i = 1 beyond it, for rocky 
and icy planetesimals respectively. 

The surface density of planetesimals evolves as a result of 
accretion and ejection by the forming planets. The same proce- 
dure as in A05 is adopted. Planetesimals that can be accreted by 
a growing planet are those within the planet's feeding zone, here 
assumed to be an annulus of 5Ru at each side of the orbit. The 
solids surface density inside the feeding zone is considered to be 
constant, i.e. the protoplanet instantaneously homogenises it as a 
result of the scattering it produces among planetesimals. Ejected 
planetesimals are considered to be lost. 

We stress that this work is the first one of a series of papers. 
Here, effects as planetesimal drifting due to gas drag, fragmen- 
tation and planetesimal size distribution are neglected, but will 
be included in future works. 



2.3. Gas accretion: attached phase 

Equations 



Planetary growth proceeds through solids and gas accretion. 
Gas accretion is a result of a planet's contraction, and is com- 
puted by solving the standard internal structure equations: 



dr 3 
d~M r 
dP 
~dM r 
dT 
dP 



3 

4np' 

-G(M r + M core ) 
47rr 4 

V a d or V rad , 



(10) 

(11) 
(12) 



where r, P, T are respectively the radius, the pressure and the 
temperature inside the envelope. These three quantities depend 
upon the gas mass, M r , included in a sphere of radius r. The 
temperature gradient is given by the adiabatic (V a d) or by the ra- 
diative gradient (V ra d), depending upon the stability of the zone 
against convection, which we check using the Schwarzschild cri- 
terion. These equations are solved using the equation of state 
(EOS) of Saumon et al. (1995). The opacity, which enters in the 
radiative gradient V ra d is computed from Bell and Lin (1994). 
In this work we assume that the grain opacity is the full in- 
terstellar opacity. However, Podolak (2003) and Movshovitz & 
Podolak (2008) showed that the grain opacity in the envelope 
of a forming planets should be much lower than the interstellar 
one. Reducing the grain opacity accelerates the formation of gi- 
ant planets (Pollack at al. 1996, Hubickyj et al. 2005) because it 
allows the runaway of gas to start at smaller core masses. Since 
the objective of the present paper is to explore the consequences 
of the planet-planetesimal interactions, we only consider a single 
opacity, corresponding to the full interstellar opacity. 

Note that we have omitted the energy equation that gives 
the luminosity, which itself enters in the radiative gradient 
(in the parts of the planet that are stable against convection), 
and in the determination of the planet stability to convective 
motions. Including the energy equation in its standard form 

-dS 



£ m - T % (the first term resulting from the accretion of 



dL 
dM, 

planetesimals, the second one to the contraction of the planet) 
brings usually numerical difficulties. Here we follow Mordasini 
et al. (2012b) to calculate the luminosity in an easier way. Note, 
however, that we improved their approach to take into account 
analytically the energy of the core. 

The total luminosity is given by L = L cont + L m , L cont being 
the contraction luminosity and L m the accretion luminosity. We 
assume L cont to be constant in the whole planetary envelope. L cont 
is computed as the result of the change in total energy of the 
planet between two time steps t and t + dt: 



Leant — 



£ tot (f + dt) - EU0 ~ E 



gas.acc 



dt 



(13) 



where E tot is the total planetary energy and £g as , a cc = dtM g:is Ui nt 
is the energy advected during gas accretion (H; nt being the inter- 
nal specific energy). £ gas ,acc is negligible compared to the other 
terms. The luminosity due to accretion of planetesimals is 



M cnre M ( 



core J " core 



(14) 



However, the energy at the time t + dt is not known before 
computing the internal structure at this given time. To circum- 
vent this problem, we use the following approach: the energy is 
split in two parts, one related to the core, one to the envelope. 

The core energy is given by E core = -(3/5)GM;? ore /,ft core , 
the core density being assumed to be uniform. The enve- 
lope energy is assumed to follow a similar functional form: 
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£env = -^env-^envg, where g is a mean gravity, taken to be 
G(M core /R COK + M tot /^pi an et). This last formula defines k em , in 
which all our ignorance of the internal structure is hidden. In 
order to calculate the envelope energy at time t + dt, the value 
of k em is first taken to be the value resulting from the structure 
at time t. Then, iteration on k em is performed until convergence 
is reached. In general, only a first order correction is enough to 
reach a satisfactory solution. 
Boundary conditions 

The internal structure equations are solved for M r varying 
between the core mass M core , and the total planetary mass. Four 
boundary conditions are given, namely the core radius, the total 
planetary radius, and r sul -f jP i anet and Psurf.pianet; the temperature 
and pressure at this point. Given the boundary conditions, the 
differential equations have only one solution for a given total 
planetary mass. 

The core radius is given as a function of the core luminosity 
and the pressure at the core surface by the following formula, 
which constitute a fit to the results of Valencia et al. (2010): 



9800km 



x 0.28+0.02 
iKi core 1 



10 



-[l°g 10 (l+wVIr))/7] : 



(15) 



where P core is the pressure a the core-envelope interface, ex- 
pressed in GPa. 

The radius of the planet, R M , is given by (Lissauer et al. 
2009): 



R 



M 



GM 



h + k 2 R H 



(16) 



where cl is the square of the sound velocity in the disc midplane 
at the planet's location, k\ — 1 and £2 = 1/4. At the planet's 
surface, the temperature and pressure are given by: 



surf.planet 



and 



1 disc + 



3tL 



\6no-pR 2 



1/4 



(17) 



(18) 



/'surf, planet — /'disc 5 

with t = A-(7di sc ,pdisc)Pdisc^Mi L I s ^ e planet luminosity, and 
Tdisc, Pdisc, /'disc are the temperature, density and pressure in the 
disc midplane at the location of the planet. 



- the pressure, that includes the dynamical pressure due to gas 
free falling from the disc to the planet, 



/'surf.planet — /'disc 



M gas 
4^ 



v«. 



(19) 



In this equation, Vg is the free falling veloc- 
ity from the Hill radius to the planetary radius, 
v ff = - a/2GM x (1/R M - 1//?h)- Note that the plane- 
tary radius is not known a priori, but computed as a result 
of integrating Eqs. ( fT0] > to ( fT2| >. 



- the maximum accretion rate, M„ 



which is equal to 



= max [F(a M + Rh) , 0] +min [F (a M - R H ) , 0] (20) 



where F = 37rvE + 6nr^ is the mass flux in the disc. 
Geometrically, the maximum accretion rate that can be pro- 
vided by the disc is equal to the mass flux entering the 
planet's gas feeding zone. The gas can enter either from the 
outer parts of the disc (which is the general case), or from 
the inner part of the disc (which can be the case in the outer 
part of the disc). 

Therefore, during a time step, a planet has access to the mass 
delivered at its location by the disc (M gas max x dt), and a mass 
reservoir made of the gas mass already present in the planet's gas 
feeding zone (see also Ida and Lin 2004). This reservoir of gas is 
assumed to be empty when the planet is massive enough to open 
a gap (which coincides with the transition to type II migration, 
see next section). However, the feeding zone continues receiving 
gas due to viscosity at the local accretion rate (Eq. [20|>. 

2.4. Orbital evolution: disc-planet interaction 

Disc-planet interaction leads to planet migration, which can oc- 
cur in different regimes. For low mass planets, not massive 
enough to open a gap in the protoplanetary disc, migration oc- 
curs in type I (Ward 1997, Tanaka et al. 2002, Paardekooper et 
al. 2010, 201 1). For higher mass planets, migration is again sub- 
divided in two modes: disc-dominated type II migration, when 
the local disc mass is larger than the planetary mass (the mi- 
gration rate is then simply given by the viscous evolution of the 
protoplanetary disc), and planet-dominated type II migration in 
the opposite case (see M09). The transition between type I and 
type II migration occurs when 



2.3.1. Gas accretion: detached phase 

By solving the differential equations ( fT~0| > to ( p"2| ) with the bound- 
ary conditions mentioned above, one can derive the planetary 
envelope mass as a function of time, and therefore the gas ac- 
cretion rate M gas . However, the rate of gas accretion that can be 
sustained by the protoplanetary disc is not arbitrary, and is in 
particular limited by the viscosity. When the gas accretion rate 
required by the forming planet is larger than the one that can be 
delivered by the disc, M gas max , the planet goes into the detached 
phase. 

In the detached phase, the planetary growth rate by gas ac- 
cretion does not depend upon its internal structure, but is rather 
given by the structure and evolution of the disc. During this 
phase, the internal structure is given by solving the same equa- 
tions ( p~0] > to ( [12} , this time for a mass M r ranging from M core to 
A^pianet (which is known). The boundary conditions are the same, 
except two of them: 



3 //disc 

4~kV 



50M* 
MRe 



1, 



(21) 



(Crida et al. 2006), where //disc is the disc scale-height at the lo- 

cation of the planet, and Re = - M — is the macroscopic Reynolds 
number at the location of the planet (v is the same as the one 
used for the disc evolution). 

First models of type I migration (Ward 1997, Tanaka et al. 
2002) predicted so rapid migration rates that it was necessary to 
reduce arbitrarily the migration rate by a constant factor, named 
fi in A05 and M09, in order to reproduce observations. Since 
these first calculations, type I migration has been studied in great 
details, and new formulations for type I migration rates are now 
available (Paardekooper et al. 2010, 201 1). We use in our model 
an analytic description of type I migration, which reproduces the 
results of Paardekooper et al. (2011). A detailed description of 
this model is presented in Dittkrist et al. (in prep.), and prelimi- 
nary results have been presented in Mordasini et al. (2010). 
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3. The accretion rate of solids 

The growth of the solid component of a protoplanet, M core , is 
assumed to be due to the accretion of planetesimals. Adopting 
the particle-in-a-box approximation, its growth rate can be cal- 
culated with 



dM cl 



df 



'27rZ m R 2 H 

K ^orbital 



colb 



(22) 



(Chambers 2006), where Y. m is the solids surface density at the 
location of the protoplanet and P or bitai is its orbital period. The 
collision rate, P co u, is the probability that a planetesimal is ac- 
creted by the protoplanet. This probability depends upon the rel- 
ative velocity between planetesimals and the protoplanet which, 
in turn, depends upon the planetesimals eccentricities and incli- 
nations. We represent by e (i) the root mean square of the eccen- 
tricity (inclination) of planetesimals. Planetesimals are found to 
be in different velocity regimes depending on their random ve- 
locities. These regimes are known as high-, medium- and low- 
velocity regime. Each regime is characterised by a range of val- 
ues of planetesimal's reduced eccentricities (e = ae/R H ) and 
inclinations (i = ai/Ru): the high-velocity regime is defined by 
e, i > 2, the medium-velocity regime by 2 > e, i > 0.2 and the 
low-velocity regime by e, i < 2. This leads to different collision 
rates 



high 



(R + r m ) 2 I 6R H I G (fi) \ 

^^ |/fW+ (/^)# 



med — 



2ttR 2 h 
(g + r m f 
4nR 2 J 



17.3 



232fl H 
R + r m 



Plow - 11-3 



R + r„ 
Rh 



1/2 



(23) 
(24) 
(25) 



(see Inaba et al. 2001 and references therein) where R is the ra- 
dius of the protoplanet (in the case of a solid body without a 
gaseous envelope R is its geometrical radius), r m is the radius of 
the planetesimals, ft = i/e, and the functions I F and I c are well 
approximated by 



1 + 0.95925/3 + 0.7725 1/3 2 
/?(0.13142 + 0.12295/3) ' 
1 + 0.3996/3 
/3(0.0369 + 0.048333/3 + 0.006874/3 2 ) ' 



(26) 
(27) 



for < P < 1, which is the range of interest for this work 
(Chambers 2006). 

According to Inaba et al. (2001) the mean collision rate can 
be approximated by 



P co „ = min (P med , (P h 2 gh + P^T m ) . 



(28) 



When an embryo is able to gravitationally bind gas from 
its surroundings it becomes more difficult to define its radius, 
which is not just the core radius. For the purpose of the colli- 
sion rate, the capture radius of the protoplanet should depend 
upon the mass of the protoplanet, upon the planetesimals' veloc- 
ity with respect to the protoplanet, upon the density profile of 
the envelope, p(r), and upon the size of the accreted planetesi- 
mals (smaller planetesimals are more affected by the gas drag of 
the envelope and therefore are easier to capture). As in Guilera 
et al. (2010), here we adopt the prescription of Inaba & Ikoma 



(2003) where the capture radius R can be obtained by solving the 
following equation 



_ 3 p(R)R 

r m — ~ 

2 p m 



v 2 el + 2GM(R)/R 
v 2 el + 2GM(R)/R H 



(29) 



where p m is the planetesimals' bulk density, G is the gravitational 
constant and the relative velocity v re i is given by 



Vrel 



= v k V5/8e 2 + 1/2 i 2 , 



(30) 



with Vk is the keplerian velocity (vt = Cla). This simple formula 
for the capture radius approximates well more complex models 
(as the one described in A05) with the advantage that it reduces 
the computational time. 

It is clear from the above equations that the accretion rate of 
solids depends upon the eccentricities and inclinations of plan- 
etesimals, which define their relative velocities with respect to 
the embryo: the higher the relative velocity is, the less likely 
planetesimals are captured by the embryo. The eccentricities 
and inclinations of planetesimals are affected by the damping 
produced by the nebular gas drag, by the gravitational stir- 
ring of the protoplanet (protoplanet-planetesimal interactions) 
and, to a lesser extent, by their mutual gravitational interactions 
(planetesimal-planetesimal interactions): 



de 2 
~d7 

df 
dt 



&<P_ 

df 
dt 



drag 



drag 



de 2 
' ~d7 

df 
dt 



VS,M 



de 2 

df 
dt 



VS,m 



(31) 
(32) 



The first term represents the effect of the nebular gas drag, the 
second term the viscous stirring produced by an embryo of mass 
M and the third term the planetesimal-planetesimal viscous stir- 
ring. 

The drag force experienced by a spherical body depends 
upon its relative velocity with respect to the gas. If we con- 
sider that the protoplanetary nebula is mainly composed by H2 
molecules, the mean free path of a molecule of gas is 



a = (HH 2 o- H2 r 



(33) 



where « H2 is the number density of H2 molecules and cr H2 is the 
collision cross-section of an H2 molecule. Depending upon the 
ratio between the planetesimal's radius and the mean free path 
of the molecules, three drag regimes can be defined (Rafikov 
2004 and references therein). The first two drag regimes are for 
planetesimals which radii are larger than the mean free path, 
r m > A. These are the quadratic and the Stokes regime. For 
the distinction of these regimes we adopt the criterion proposed 
by Rafikov (2004) in terms of the molecular Reynolds num- 
ber Re mo ] = v re ir m /v mo i, where v mo i is the molecular viscosity, 
v m o\ = Ac s /3. If Re mo i > 20 we assume that the gas drag is in the 
quadratic regime, and the differential equations for the evolution 
of the eccentricity and inclination are given by 



de 2 
df 

df 
dt 



drag 



drag 



2e 2 19 

Tdrag \4" 47T 



J + ^2 + if 
71 



1/2 



T drag 



yf + L e 2 + if 



1/2 



(34) 
(35) 
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(Adachi et al. 1976 corr ected by Inaba et al. 2001 fj with 
£ a 1.211. The value of 77 depends upon the distance to the star, 
on the gas density and on the pressure gradient, dP/da, 

( ) " 1 dP (36) 
2£2 2 ap gas da ' 

where p eas and dP/da are derived from the disc model. 
The gas drag timescale is 

8pm Fm 



^"drag — 



(37) 



3C D p g a S Vk 

where Cd is the drag coefficient, which is of the order of unity. 

The Stokes regime occurs for r m > A and Re mo i < 20, and the 
equations for the eccentricity and inclinations of planetesimals 



(38) 



are 






de 2 




3 Ac s p gas e 2 


dt 


drag 


2 P m r 2 m 


di 2 




3 Ac s p gas f 


dt 


drag 


4 Pm r l 


(Adachi et al. 1976, Rafikov 2004) 



(39) 



Note that in this paper, as for example in Stepinski & 
Valageas (1996), we defined two different Reynolds numbers, Re 
and Re mo i, and two different viscosities, v and v mo i. The macro- 
scopic quantities (Re and v) are a measure of the fluid dynamics 
of the disc in a global scale (to compute the evolution of the disc) 
while the microscopic quantities (Re mo i and v mo i) characterise 
the local state of the gas and are used to calculate the eccentrici- 
ties and inclinations of planetesimals. 

When r m < A, the third regime, Epstein regime, takes place, 
and eccentricities' and inclinations' evolution follow the equa- 
tions 

^ = - e 2 ^i, (40) 

= 4—- w> 

drag 1 Pm r m 



dt 

df 

dt 



drag 



(Adachi et al. 1976, Rafikov 2004). In this paper we consider 
that the population of planetesimals is represented by spherical 
bodies of a single size. It is worth mentioning that, although 
we allow for the three drag regimes according to the above- 
mentioned criterions, for the ranges of planetesimal sizes con- 
sidered in this work (100-0.1 km) and for the kind of interaction 
we are mostly interested in here (protoplanet-planetesimal inter- 
actions), planetesimals are found to be mainly in the quadratic 
regime. Therefore, in most of the cases, for determining the 
solids accretion rate the effect of the gas drag is governed by 
Eqs. ((34) - ((35). 

Planetesimals' eccentricities and inclinations are excited by 
the presence of a protoplanet. Ohtsuki et al. (2002) studied the 
evolution of the mean square orbital eccentricities and inclina- 
tions and introduced semi-analytical formulae to describe the 
stirring produced by the protoplanet 



de 2 
~dt 

df 
dt 



VS.M 



VS.M 



M 



3/>M*P orbital 
M 

3/W*P orbital 



'vs, 



Q 



vs> 



(42) 



(43) 



3 We did not find any noticeable difference in our results when using 
the original formulas of Adachi et al. (1976) or the corrected ones by 
Inaba et al. (2001). 



where b is the full width of the feeding zone of the protoplanet 
in terms of their Hill radii (here we adopt b ~ 10), and Pys and 
2vs are given by 

ln(l + 10A 2 /e 2 ) + 



10A 2 
72/pvsOS) 



Q 



vs 



4l 2 + 0.2le 3 



+ 



10A 2 e 
72/ Q vsOS) 



ln(l + A ), 



ln(l + 10A 2 e 2 ) + 



(44) 



ln(l + A 2 ), (45) 

nei 

with A = l(e 2 + P)/12. The functions 7 PV s(/3) and 7qvs(/3) can be 
approximated, for < /3 < 1 , by 

/3- 0.36251 



/pvsOS) - 
/qvsOS) - 



0.061547 -1- 0.16112/3 + 0.054473/3 2 ' 
0.71946 -/3 



(46) 
(47) 



0.21239 + 0.49764/? + 0.14369/3 2 ' 

(Chambers 2006). The excitation that the protoplanet produces 
on the planetesimals weakens with the increase in the distance 
between the protoplanet and the planetesimals, i.e. further away 
planetesimals are less excited. Here we follow the approach of 
Guilera et al. (2010) and consider that the effective stirring is 
given by 



(48) 



de 2 


eff 


de 2 




~dJ 








VS,M 


VS.M 


df 


eff 


d; 2 




dt 




= ' (A) d7 




VS.M 


VS,M 



(49) 



where /(A) ensures that the perturbation of the protoplanet is 
confined to its neighbourhood, 

-1 



/(A) 



1 + 



A 
riRn 



(50) 



with A = \om - a m \, where om is the semi-major axis of the 
protoplanet and a„, is the semi-major axis of the planetesimal. 
Although the functional form is arbitrary, the scale on which 
the stirring acts is similar to the one found in A^-body calcu- 
lations (excluding the effects of resonances). For this work we 
have chosen n — 5 to limit the perturbation of the planet to its 
feeding zone. In the future, with the aid of A^-body calculations, 
we plan to have a better semi-analytical function to characterise 
the extent of the planetary perturbation. 

We also consider that planetesimals' eccentricities and incli- 
nations are stirred by their mutual interactions. For a population 
of planetesimals of equal mass m, the evolution of their eccen- 
tricities and inclinations are well described by 



de 2 
df 

df 
dt 



VS.m 




VS,m ' * 

(Ohtsuki et al. 2002), with 

v 1/3 

/ 'm 

h„ 



VS, 



3M* 



(51) 
(52) 

(53) 
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In this case Pys and (2vs are evaluated with the reduced ec- 
centricity and inclination relative to the planetesimal mass (i.e. 
e — 2e/h m , i = 2i/h,„). Note that there is no dynamical fric- 
tion term in Eqs. ( |5T| - ( |52| , as it vanishes when a single mass 
population of planetesimals is considered (Ohtsuki et al. 2002). 
Although, strictly speaking we have two populations of planetes- 
imals (rocky and icy bodies depending if they are inside or be- 
yond the ice line) that have the same size but not the same mass 
(because of the difference in their density), the region where the 
two types of planetesimals are present at the same time is very 
narrow. 

In this work we neglect changes in the mass of planetesimals 
due to fragmentation and changes in the surface density owe to 
planetesimal drifting. 

3.1. Comparison with the previous solids accretion rate 

In previous works (A05, M09), the solids accretion rate has been 
treated in a very simple way leading to an underestimation of the 
formation timescale of planets. In those works, the prescription 
for planetesimals' eccentricities and inclinations was the same 
as in Pollack et al. (1996), were it was assumed that planetesi- 
mals' inclinations depend only on planetesimal-planetesimal in- 
teractions. Under this assumption, the reduced value of planetes- 
imals' inclination, i, was prescribed as: 



/ = 



(54) 



where ve is the escape velocity from the surface of a planetes- 
imal. This means that planetesimal' s inclination i - i Rn/a is 
constant independently of the mass of the planet. On the other 
hand, eccentricities are assumed to be controlled by both plan- 
etesimals and protoplanet stirring, its value given by: 



e = max(2i, 2). 



(55) 



Therefore, if e = 2i the protoplanet is growing according to the 
runaway regime because e and i would be independent of the 
mass of the protoplanet. If e = 2, the eccentricity of planetes- 
imals would be affected by the presence of the protoplanet, so 
to a certain extent the stirring of the embryo is taken into ac- 
count. However, this condition corresponds to the protoplanet- 
planetesimal scattering in the shear-dominated regime. Ida & 
Makino (1993) showed that the shear-dominated regime lasts 
for only a few thousand years, after which planetesimals are 
strongly stirred by the protoplanet. During the shear-dominated 
period, eccentricities and inclinations of planetesimals in the 
vicinity of the protoplanet remain low. This leads to an accre- 
tion scenario which is much faster than that corresponding to the 
oligarchic regime (the oligarchic regime usually occurs in the 
dispersion-dominated regime). To show clearly the difference 
between this, let's call it "quasi-runaway" accretion of solids, 
and the oligarchic regime, we performed two simulations that 
are identical in all parameters except for the prescription of e 
and i. The planet is assumed to form in situ at 6 AU. Accreted 
planetesimals are 100 km in radius. No disc evolution is con- 
sidered, so simulations are stopped when the planet reaches one 
Jupiter mass. For the quasi-runaway regime we use Eqs. ( [55} and 
p4| i to calculate e and i, while for the oligarchic growth we solve 
the differential equations presented in the previous section. Fig. 
[2] (left panel) shows the ratio of oligarchic and quasi-runaway 
eccentricities and inclinations as a function of the mass of the 
planet. In the case of the eccentricity, the one corresponding to 



the oligarchic regime is ~ 4 times larger than in the quasi run- 
away regime. The oligarchic inclination is several tens of times 
bigger than the quasi-runaway one. As a consequence, the ac- 
cretion rate of solids is much smaller in the oligarchic than in 
the quasi runaway regime because in the oligarchic regime plan- 
etesimals are more exited and are more difficult to accrete. A 
comparison between the two accretion rates of solids is shown 
in Fig.[2](right panel). Due to the smaller accretion rate in the oli- 
garchic regime, while in the quasi runaway regime it takes less 
than 1 Myr to form the planet, in the oligarchic regime formation 
is much longer, taking 3.25 x 10 7 years. 

4. Results 

In the previous section we introduced the main characteristics of 
the planet formation model, with special focus on the differences 
in the physical and numerical model with respect to A05. In the 
following sections we will concentrate on the impact that the ac- 
cretion rate of solids has on the formation of giant planets. Here 
we will consider the formation of isolated planets, i.e. only one 
planet per disc. The computations of planetary system formation 
will be presented in other papers (Alibert et al. 2012, Carron et 
al. 2012). 

As we described in Sect. [3] the treatment of the evolution of 
eccentricities and inclinations of planetesimals intends to min- 
imise the assumptions on their values (keeping in mind that it is 
not an A^-body calculation, but the adopted formulas reproduce 
iV-body results of planetesimals accretion rates and excitation). 
In order to consider a realistic accretion rate but not too com- 
putationally expensive, Thommes et al. (2003) considered that 
planetesimals' eccentricities and inclinations can be estimated 
assuming that the stirring produced by the protoplanet is instan- 
taneously balanced by the gas drag. The approximation to the 
equilibrium values of e (e e q) can be derived by equating the stir- 
ring timescale and the damping timescale, resulting in 



1.7 



OT l/lS M l/3 p 2/15 



(56) 



The equilibrium value for i (/ eq ) is assumed to be half the 
value of e eq , as this relationship has been shown to be a good 
approximation in the high-velocity cases (Ohtsuki et al. 2002): 



1 



(57) 



However, it is not clear whether planetesimals are always in 
equilibrium, especially if we consider that depending on their 
mass, planetesimals are differently affected by gas drag, and that 
during its formation, a planet migrates and the protoplanetary 
disc evolves. On the other hand, if equilibrium is attained, it is 
interesting to compare the equilibrium values obtained by solv- 
ing explicitly Eqs. (3T)-([32| with the approximations given by 
Eqs. (|56|>-(|57)i based on timescales considerations. 

Fortier at al. (2007, 2009) and Benvenuto et al. (2009) 
assume, for simplicity, the equilibrium approximation of 
Thommes et al. (2003) in their in situ, giant planet formation 
models. However, Chambers (2006) finds that important devia- 
tions from equilibrium occur at the very beginning of the growth 
of the embryo but eventually equilibrium is attained. In the cases 
he shows, these deviations do not seem to have a noticeable ef- 
fect in the final mass of the planet as long as no time restriction 
for the lifetime of the disc is assumed. On the other hand, Guilera 
et al. (2010, 201 1) do not use any approximations and calculate 
explicitly e and i by solving the corresponding time evolution 
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Fig. 2. Left panel shows the ratio between eccentricities (red) and inclinations (blue) as a function of the planet's mass considering 
both the oligarchic and the quasi runaway growth. Clearly, planetesimals excitation is several times higher in the oligarchic growth 
than in the quasi-runaway growth. Right panel shows the impact of planetesimals excitation in the accretion rate of solids: in the 
case of the quasi-runaway growth (yellow) the accretion rate is much larger than in the oligarchic growth (green). 



differential equations in their giant planet formation model. No 
study comparing the equilibrium approximation and the explicit 
calculation of e and i has been made using a self consistent gi- 
ant planet formation model. Moreover, out of equilibrium effects 
can be important, not only at the beginning of the formation of a 
planet. When planets migrate they can enter regions where plan- 
etesimals are in principle cold (low values of e and /), or already 
excited if there is another planet growing in the neighbourhood. 
Depending on the ratio between the stirring and the migration 
timescale, one can expect some cases where the equilibrium ap- 
proximation may not be accurate. 



We study the formation of giant planets considering a self- 
consistent model for the interplay between the disc evolution, 
accretion by the growing planets and gas driven migration. In 
this paper, we focus our study on the formation of single planets, 
looking in detail at the dependence of planetary growth upon the 
planetesimal size and the differences in the final results between 
the equilibrium approximation and the explicit calculation of e 
and i, all together with planetary migration. We do this analysis 
in four steps. First we consider the formation of 1 M ffi planet, 
neglecting the presence of an envelope, of planetary migration 
and of disc evolution (Sect. |4.1| ). Second, we compute the full 
formation of a planet (except migration), which is assumed to 
be over when the gaseous component of the disc disappears. We 
keep the in situ formation hypothesis to ease the analysis (Sect. 
|4.2| >. Third, we allow for gas driven migration during the forma- 
tion of the planet (Sect. |4~3] >. Fourth, we generalise the examples 
presented in the third step by considering a wide range of plau- 
sible protoplanetary discs and initial locations for the embryo to 
have an overview of all possible outcomes (Sect.|4~4|. 



4.1. Formation of 1 M e planet 

We analyse first the formation of 1 M e planet (we stop the cal- 
culation when this mass is reached), neglecting the presence of 
an envelope, planetary migration and disc evolution. The point of 
this case is to focus on the initial stages of the accretion and anal- 
yse the importance of the size of the accreted planetesimals and 
the implications for the growth of the embryo when considering 
the equilibrium approximation. It is because we want to show 
clearly the consequences of these assumptions that we neglect 
other physical processes that act simultaneously. This means, 
for example, that for a fixed mass of the embryo and indepen- 
dently of the elapsed time, the state of the protoplanetary disc is 
the same in terms of surface density (solids and gas). The same 
applies for the capture radius of the planet. In this example, the 
embryo is assumed to be located at 6 AU, where the initial solids 
surface density is E,„ = 10 g cirT 2 and the density of the nebular 
gas is p gas = 2.4 x 10~ 9 gcirT 3 . For this disc the snow line is at 
3.5 AU. The initial surface density profile of the disc is given by 
Eq. ([H} with y = 0.9 and ac = 127 AU. The initial mass of the 
embryo is 0.01 M®. For the equilibrium approximation we adopt 
Eqs. (56} - ( [57} to calculate the values of e and i. For the explicit 
calcu ation of the eccentricity and inclination of planetesimals 
we solve Eqs. pi) - ( [32| . To solve Eqs. ( |3"T) - ( [32] i, initial condi- 
tions for e and ; must be given. We consider two possibilities for 
the initial conditions that we think bracket the parameter space. 
On one hand, we consider that the planetesimal disc is initially 
cold, and planetesimals' eccentricities and inclinations are given 
by the equilibrium value between their mutual stirring and the 
gas drag. These values can be derived by equating the stirring 
timescale and the damping timescale, that results in 



= 2.31 



m 4/15 2 l/5 a l/5 p . 



2/15 



(58) 
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Fig. 3. In situ formation of 1 M e planet at 6 AU for different radii of accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; 
orange, 0.1 km). Solid line is for the equilibrium approximation of eccentricities and inclinations (see Eqs. p6| ) - (|57j)). Dashed and 
dotted lines are for the explicit calculation of e and i solving the differential equations Eqs. pTh - (32 1. Two different initial values, 
<?o, z'o> were considered: dashed line represents the case where <?o and j'o are given by Eqs. (|59]l - (|58), while dotted line assumes <?o, *o 
to be the initial equilibrium values of the equilibrium case. Top-left panel depicts the mass growth of the protoplanet as a function of 
time. Top-right panel shows the time evolution of the eccentricity. Bottom-left panel shows the eccentricity as a function of the mass 
of the embryo. Clearly, in all cases, after an out of equilibrium state, equilibrium is attained. Bottom-right panel plots the accretion 
rate of solids as a function of the mass of the embryo. Note that, as a function of the planet mass, the difference in the accretion rate 
is only a consequence of the difference in e and i. 



•m-m _ m-m 

'eq — 2 e 1 



(59) 



This means that we are assuming that the embryo instanta- 
neously appears in an unperturbed planetesimal disc. The other 
extreme situation is to assume a hot disc, where planetesimals 
are already excited by the embryo and their initial eccentricities 
and inclinations are those corresponding to the value of equilib- 
rium between the stirring of the embryo (0.01 M ffl ) and the gas 
drag, approximated by Eqs. ( |56| ) - (|57) . The initial values of e 
for these cases are given in table |2| where e™~"' correspond to 
Eq. (|58ja 
to be e/2 



eq o to Eq. ( 56 1. The initial values of i are assumed 
= eo/2). 



We perform calculations for the equilib- 
rium approximation and the two sets of initial conditions for the 
explicit calculation, for four radii of the accreted planetesimals: 
100, 10, 1 and 0.1 km. 

Fig.|3]shows the results of the simulations mentioned above. 
The top-left panel depicts the mass growth of the planets as a 
function of time. The first evident thing from this plot is the 
timescale difference in the formation of 1 M e planet depend- 



Table 2. Two sets of initial values of the eccentricity for the ex- 
plicit calculation of its evolution. 



r m 


gfn—m 
cq.O 


e eq,0 


100 km 


8.4 xl0~ J 


2.13 xlO- 2 


10 km 


1.3 xl0~ 3 


1.35 xl0~ 2 


1 km 


2.0 xl0~ 4 


8.70 xl0~ 3 


0.1 km 


3.3 xlO" 5 


5.53 xl0~ 3 



ing on the size of the accreted planetesimals: while for accreted 
planetesimals of 100 km it takes ~ 10 7 yr to grow from 0.01 
M ffl mass to 1 M ffl , it takes ~ 10 5 yr in the case of 0.1 km plan- 
etesimals. This difference in the growth timescale is entirely due 
to the fact that large planetesimals are less damped by gas drag 
than the smaller ones. Hence, while they are stirred up by the 
massive embryos to about the same value, they keep larger e and 
i values (bottom-left panel) making the accretion process much 
slower. Keep in mind that the disc does not evolve and planets 
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do not migrate, so differences depend only on the planetesimal 
size. Note that for a fixed mass of the embryo the accretion rate 
of solids differs in two orders of magnitude between the two ex- 
treme cases (bottom-right panel). 

If we now turn our attention to the differences in growth 
rate for a fixed planetesimal size, but for different approaches in 
the calculation of the eccentricities and inclinations, we see that 
adopting the equilibrium approximation may lead to significant 
differences in the mass of the planet, more evident when the ac- 
creted planetesimals are small (top-left panel, compare solid line 
with dashed or dotted lines of the same colour). For 100 km and 
10 km planetesimals we do not see deviations from equilibrium 
when the initial conditions for e and i are given by Eqs. (|56i 



- ( |57| (red and blue dotted lines, top-right panel). If the initial 
values for e and i are given by Eqs. (58 i - ( |59] > the equilibrium 
values given by Eqs. (56 1 - (57 1 are reached in an almost negli- 
gible fraction of the formation timescale. We conclude that the 
equilibrium approximation for e and i gives results that nicely 
agree with more complex calculations, as far as the planetesimal 
size is relatively big or the evolution time is sufficiently long. 
However, as shown on the same picture, the growth of a small 
planet is very long when considering such massive planetesi- 
mals, and the formation of a gas giant under these conditions is 
highly compromised. 

The situation becomes more critical when we consider 
smaller planetesimals. Deviations from equilibrium are evident 
regardless the initial values adopted for e and i, especially for 
r m = 0.1 km. Nevertheless, equilibrium is always attained, al- 
though the equilibrium values for e and i are lower than those 
given by Eqs. po*) - ( |57] i, especially in the case of the inclination. 
As we can see from Fig. [4] the approximation yS eq = i/e = 1/2 is 
not very good for smaller planetesimals. Note that the real equi- 
librium value of f3 depends upon the planetesimal size. While 
jS eq = 1 /2 is a good approximation for planetesimals larger than 
~ 1 km, it is not the case for smaller planetesimals, which tend to 
have lower inclination values than 1/2 e. We find that for small 
planetesimals, the velocity regime is in the limit between the 
high and medium velocity regime (e , i < 2), therefore eccentric- 
ities are more effectively excited than inclinations (see Ohtsuki 
et al. 2002). The low values of <?, i and ft increase the accretion 
rate, as can be seen in the bottom-right panel of Fig. [3] speed- 
ing up the formation of the embryo (relative to the equilibrium 
approximation case) . 

4.2. Planet formation without migration 

In this section we analyse the complete formation of a planet for 
the same cases as before. Now, the planet accretes mass (solids 
and gas) as long as the disc is still present. This means that we 
allow for disc evolution (Sect. 2.1 1, therefore planet formation is 
considered to be over when the disc disperses. The presence of 
the gaseous envelope is taken into account for the calculation of 
the capture radius. However, we still neglect the migration of the 
embryo to ease the analysis of the dependence on planetesimals' 
e and i. 

As it has been already shown in other works, the presence of 
the envelope increases the capture radius, speeding up the for- 
mation of the planet. The enhancement in the capture cross sec- 
tion makes the accretion of solids more effective. This effect is 
more noticeable for small planetesimals. Compared to the cases 
of Sect. ( |4.1| i, the formation timescale of a 1 M e can be reduced 
up to ~ 35 % (for the smallest planetesimals, r m =0.1 km), due 
only to the enhancement in the capture cross section. The pres- 
ence of an atmosphere, even if its mass is negligible compared 
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Fig. 4. Inclination-eccentricity ratio for the out of equilibrium 
cases of Fig. [3] Red is for 100 km planetesimals, blue for 10 km, 
green for 1 km and orange for 0. 1 km. Dotted line represents the 
case of an initially hot disc while dashed line of an initially cold 
one. The grey line is the standard equilibrium value, /3 eq = i/e = 
1/2. 



to the total mass of the planet, has to be considered for embryos 
as small as 0. 1 M ffl as it plays an important role for the accretion 
of solids. 

Fig. 5] is the analog of Fig. [3] for the complete formation 
of the planets. As the disc disperses after 6 Myr, in the cases 
where the accretion of planetesimals is slow (r m = 100, 10 km), 
the final masses of the planets are lower than 1 M e . The differ- 
ences in growth observed in the previous section as a function 
of the planetesimal size have dramatical consequences for the fi- 
nal mass of the planet, which can be ~ 0.1 M e if the accreted 
planetesimals have a radius of 100 km, ~ 0.8 M e for planetes- 
imals of 10 km, ~ 1200 M e (3.7 Jupiter masses, Mj) for 1 km 
planetesimals, and 7100 M e (22 Mj) for 0.1 km planetesimals. 
These numbers evidence how differences in the accretion rate 
of solids (here regulated by the size of the accreted planetesi- 
mals, bottom-right panel) impact on the final mass of a planet 
and the non-linear aspect of planet formation in the core accre- 
tion model: once the critical mass is attained, the very rapid ac- 
cretion of gas leads rapidly to massive planets. The fact that big- 
ger planets form when the accreted planetesimals are small is a 
consequence of the gas drag, which operates in two ways that 
combine positively: nebular gas drag is more effective in damp- 
ing planetesimal's eccentricities and inclinations when planetes- 
imals are small (bottom-left panel) and atmospheric gas drag is 
able to deflect more distant planetesimals' trajectories therefore 
enlarging the capture radius of the planet. In fact, for the cases 
of 1 km and 0.1 km accreted planetesimals, embryos grow to 
become big giants. When a massive solid embryo is formed it 
triggers the accretion of gas, leading to the formation of a gas 
giant planet (green and orange lines). 

Indeed, planets can end up very massive if they enter the 
runaway phase of gas. As explained before, during the attached 
phase, the accretion of gas is a result of solving the equations 
presented in Sect. [23] The planet will remain attached to the disc 
until its accretion rate exceeds the maximum amount of gas that 
the disc can deliver. When this condition is not any more satis- 
fied, it goes into the detached phase. During the detached phase, 
the maximum accretion rate is given by Eq. ( 20 1. As in M09 and 
Ml 2, here we follow the results of Kley & Dirksen (2006) for 
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Fig. 5. In situ planetary formation at 6 AU for different radii of the accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; 
orange, 0.1 km). The same line code as in Fig.[3]is adopted. The formation of the planet ends when the gas component of the disc 
is dispersed. 



the accretion of gas during the detached phase. Planet-disc inter- 
actions can lead to eccentric instabilities, which means that the 
planet can enter regions that are outside its gap and have full ac- 
cess to the gas present in the disc, with no limitation for accretion 
except for the ability of the disc to supply it. Although it is not 
clear whether all planets can suffer from an eccentric instability, 
for the sake of simplicity we assume this is the general situation 
in our simulations. On the other hand, other works (e.g. Lissauer 
at al. 2009) include a limitation for the accretion of gas when the 
planet opens a gap in the disc, as they consider that planets are 
in circular orbits. Comparatively, this assumption leads to less 
massive planets. 

The equilibrium approximation and the explicit calculation 
of e and i also bring differences in the final mass of the planet. 
Whether for large planetesimals (100 and 10 km) the equilibrium 
approximation works fine as showed in the previous section, for 
smaller planetesimals it might not be the case. For r m = 1 km, 
the final mass of the planet in the equilibrium approximation is 
537 M®, whether with the explicit calculation considering an ini- 
tial hot disc is 950 M ffl and with an initial cold disc, 1200 M ffl . 
For r m = 0.1 km, the final mass of the planet in the equilibrium 
approximation is 4745 M®, whether with the explicit calculation 
considering an initial hot disc 6745 M ffi and with an initial cold 
disc, 7100 M®. Differences in mass are the consequence of a tim- 
ing effect. When the planet grows faster, the damping produced 



by the nebular gas is more effective because the gas density in 
a younger disc is larger than in an older one. Then, the faster 
an embryo grows, the more it profits from the nebular gas drag. 
Moreover, the cross-over mass (the mass of the core for which 
the mass of the envelope equals the mass of the core), which 
in this case is ~ 30 M e , is achieved earlier when the accretion 
rate of solids is bigger. When a protoplanet reaches the cross- 
over mass its growth is dominated by the accretion of gas, that 
is already in the runaway regime. If the protoplanet starts the 
accretion of gas when the disc is younger and more massive, it 
provides a larger reservoir of gas which translates, in the end, in 
a bigger planet. 

It is important to remark that Eqs. p6*| ) - ( |57| ) are obtained 
assuming that the stirring timescale is equal to the damping 
timescale (Thommes et al. 2003), and the equilibrium approx- 
imation (z = 0.5 e). On the other hand, Eqs. ( fJT) - ( [32] ) explic- 
itly solve the coupled evolution of e and i given a set of initial 
conditions. This means that out of equilibrium states are allowed 
(e.g. as seen at the beginning of the calculations) and equilibrium 
states are reached naturally and without a fixed ratio between e 
and i. It is clear that the explicit resolution of the differential 
equations is a better physical approach. We have included here 
calculations with the equilibrium approximation just for com- 
parison. Results show that the equilibrium approximation works 
fine for larger planetesimals, but overestimates the excitation of 
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Fig. 6. Full planetary formation, now allowing for the migration of the embryos, for different radii of the accreted planetesimals 
(red, 100 km; blue, 10 km; green, 1 km; orange, 0. 1 km). The same line code as in Fig.[3]is adopted. Embryos are initially located at 
6 AU. Left panel shows the cumulative mass of the embryo (solids plus gas) as a function of time (in logarithmic scale). Right panel 
depicts the migration path of the protoplanet (in linear scale). Note that for the accretion of 1 km planetesimals, the embryos are lost 
in the Sun (indicated with crosses in the plot). The same happens for the planet that grows by accretion of 0.1 km planetesimals, in 
the equilibrium approximation (orange solid line). When planets are lost, calculations are stopped. 



smaller ones making their accretion less effective. This delays 
the whole process of planet formation. As the planet is embed- 
ded in a disc with a finite (and short) lifetime, this delay im- 
pacts in the final mass of the planet. When solving explicitly 
the equations of e and i, initial conditions for these quantities 
are needed. This is a problem because we can not be certain 
about the state of the disc at the beginning of our calculations. 
So assumptions for the initial values can not be avoided. The 
"initially cold disc" favours the growth of a planet allowing for 
high accretion rates in the first thousands of years (Figj5| bottom- 
right panel). In the "initially hot disc" this effect is not present 
as departures from equilibrium are smaller and therefore equi- 
librium is attained faster (keep in mind that the equilibrium at- 
tained when solving the differential equations can be different 
to the equilibrium approximation for the case of small planetes- 
imals). However, final results do not strongly depend upon the 
initial choices. The differences in the final mass are of 20% and 
5% for r m = 1 km and r,„ = 0.1 km respectively. Given all the 
uncertainties of the model these differences are acceptable. 



4.3. Planet formation with migration 

We now complete our analysis including the migration of the 
protoplanets. The migration model is the one presented in 
Dittkrist et al. (in prep.) and Mordasini et al. 2010. Fig.|6]shows 
the total mass of the planet and its semi-major axis as a function 
of time for the same cases as in the previous sections. Clearly, 
the situation is very different from the in situ hypothesis. To start 
with, in all the cases the final location of the planet is far from 
its initial location. Actually, the three calculations for 1 km plan- 
etesimals result in lost planets (planets cross the inner edge of the 



disc and are considered to fall into the central star). The same is 
the case for the equilibrium approximation of 0. 1 km planetesi- 
mals. 

For the cases where planets are not lost into the star, we find 
that the final masses for accreted planetesimals of r m = 100 km 
and 10 km are independent of the way e and i were computed 
(~ 0.2 M e and 0.7 M e respectively). The final locations are also 
similar for the different considerations for e and i (~ 1.7 AU for 
both r m = 100 km and 10 km). For a fixed planetesimal size, the 
migration paths for the different assumptions on the eccentrici- 
ties and inclinations are somehow shifted in time but are similar 
if we consider them as a function of mass (Fig. [7}. 

In the case of r m - 1 km, protoplanets are lost in the central 
star when their mass is ~ 20 M ffl , most of which is in the solid 
core. At ~ 8 M @ the protoplanet undergoes inward migration 
in the adiabatic saturated regime^](see Paardekooper et al. 2010, 
201 1, Mordasini et al. 2010, Dittkrist et al.). This timescale turns 
out to be shorter than the accretion timescale. The protoplanet 
covers Aa ^ 5 AU in ~ 7 x 10 5 yr. In this time it doubles its 
mass which, however, is not big enough for the planet to enter 
type II migration. Planet migration in general slows down when 
the protoplanet is able to open a gap in the disc (type II migra- 
tion) which, as a rule of thumb, happens when the planet mass is 
~ 100 M ffl . This situation is not reached in this case, where ac- 



The regime is called "adiabatic" when the local cooling timescale 
in the disc is larger than the time it takes for a parcel of gas to make a U- 
turn on the horse-shoe orbit close to the planet in the corotation region. 
A regime is called "saturated" when the contribution of the corotation 
region to the angular momentum exchange is reduced by the the ratio 
of viscous timescale and libration timescale. In this case, the migration 
behaviour will get dominated by the angular momentum exchange at 
the Lindblad resonances. 
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cretion is too slow compared to migration, resulting in the loss 
of the forming planet. For r m = 0. 1 km the differences between 
the equilibrium approximation and the explicit calculation of e 
and i are more dramatic: adopting the equilibrium approximation 
leads to the loss of the planet in the central star (for the same rea- 
son as in the previous case, the growth rate is very slow) whether 
for the explicit calculation of e and i, although the planet ends in 
an orbit very close to the central star (~ 0.2 AU), the previous 
growth of the embryo is fast enough to enable a large accretion 
of gas. In this case the planet reaches a mass that enables it to 
switch to type II migration. The planet decelerates its migration 
speed until it stops. The final masses are: 13 Mj for in initially 
hot disc and 15 Mj for an initially cold disc. The fate of an em- 
bryo (becoming a big or a low mass planet, surviving or being 
lost in the central star), as it is shown here, depends upon the size 
of the accreted planetesimals and on the assumptions we adopt 
to describe their dynamics, as these strongly impact on the accre- 
tion rate of solids and therefore on the whole formation process 
through the regulation of the growth timescale. Here we have 
shown the interplay between the evolution of the protoplanetary 
disc, the growth of the protoplanet and the operation of migra- 
tion. Note that in all the cases we are considering the same disc, 
that globally evolves with time in the same way. Differences in 
the local evolution of the disc arise due to accretion (solids and 
gas accreted by the planet are removed from the disc) and ejec- 
tion of planetesimals (when the planet is massive enough). The 
fact that in Fig. [7J for the same mass of a protoplanet, the loca- 
tion in the disc can be different is a consequence of this inter- 
play: protoplanets reach a certain location earlier or later in the 
evolution of the disc, depending on their growth rate. The state 
of the disc and the mass of the protoplanet at that moment de- 
termines its migration rate. So, independently of what regulates 
the growth of the planets, the examples here show that planet's 
growth and migration rate are tightly coupled. 

Fig. [8] shows a comparison between the migration 
timescale (r m i e = a/\a\) and the protoplanet's growth timescale 
(i" g rowth = M/\M\) for the four sizes of the accreted planetesimals 
we have considered, solving explicitly the equations for e and 
i and adopting an initially cold planetesimal disc. For accreted 
planetesimals of 100 km and 10 km, the jump in the growth 
rate of the planets at M > 0.1 M e corresponds to the cross- 
ing of the ice line. Planetesimals in the inner region (inside the 
ice line) are denser and the gas drag is less effective on them 
(see Eqs. ([34j - ( |35] l and p7)i), therefore their random velocities 
are higher and accretion rates are lower. Also, the solids surface 
density is lower, which contributes to decrease the accretion rate. 
The peaks in the migration timescale correspond to changes in 
the sense of migratioij^] 

When the accreted planetesimals have a radius of 1 km, the 
migration timescale becomes lower than the growth timescale 
when the mass of the protoplanet is ~ 10M ffi . Planetesimals rel- 
ative velocities are very high in the neighbourhood of the pro- 
toplanet and accretion becomes more difficult as the protoplanet 
increases its mass. The protoplanet grows slowly and migrates 
fast. This situation is never reverted and the protoplanet is lost in 
the central star. Being the migration rate high and the accretion 
rate not enough to counteract it, the protoplanet migrates with- 
out limit, almost at a constant mass, until it reaches the central 
star. 



5 Depending on the regime - isothermal versus adiabatic, and satu- 
rated versus unsaturated, the migration can be inward or outward, see 
Paardekooper et al. 2010, 2011, Mordasini et al. 2010, Dittkrist et al., 
in prep. 
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Fig. 7. Mass versus semi-major axis for the cases shown in Fig. 
[6] Note that the accretion rate of solids plays a major role not 
only in the growth of the planet but also in its migration path. 



In the case of accreted planetesimals of 0.1 km, the pro- 
toplanet's migration timescale is shorter than the protoplanet's 
accretion timescale (just as in the former case) for the "criti- 
cal" mass interval of a few tens to around one hundred of Earth 
masses. However, in this case accretion proceeds fast enough for 
the protoplanet to start runaway accretion of gas before reaching 
the inner edge of the disc. Therefore the planet is able to ac- 
crete mass very fast (gas accretion dominates the growth of the 
planet), which means that it becomes massive enough to enter 
type II migration. The migration rate (in type II with gap open- 
ing) decreases as the protoplanet grows in mass. Therefore, be- 
ing in the runaway phase is what saves the planet from falling 
into the star. 

The situation can be summarised as follows: when the 
growth of the planet is dominated by the accretion of solids 
in the oligarchic regime (before gas runaway accretion), the 
growth timescale is proportional to M 1 ^ 3 (Ida & Makino 1993). 
However, for planets massive enough to trigger rapid gas accre- 
tion, the growth timescale is much shorter. On the other hand, 
migration typical time scales with M~' in type I, and is indepen- 
dent of the planet's mass in type II, being generally much longer 
than in type I. As a consequence, if a planet migrates in type I 
and is dominated by the oligarchic growth, it is very likely to 
be lost in the star. It is only if the planet succeeds to become 
massive enough to start the runaway gas accretion, and quasi- 
simultaneously enter type II migration, that it can brake before 
reaching the central star. Therefore, there is a critical mass range 
between ~ 10 — 100 M ffl : a planet in this mass range is likely to 
be undergoing inward type I migration and a decreasing growth 
rate, which in this stage is dominated by the accretion of solids. 
Indeed, being massive, the protoplanet excites the random ve- 
locities of planetesimals making its accretion difficult. Therefore 
the growth of the planet is very slow at this stage, as planetesi- 
mals can not be effectively accreted. Although this is the path- 
way to the runaway accretion of gas, the transition between be- 
ing dominated by the accretion of solids and the accretion of gas, 
even fast, is not immediate. On the other hand, inward type I mi- 
gration in this mass range is still very fast. Therefore, if planets 
enter the runaway gas phase they would probably grow massive, 
change to type II migration and avoid a fatal end in the central 
star. If they remain very small, type I migration is not very harm- 
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Fig. 8. Protoplanet's growth and migration timescales (T grow th = M/\M\ and r m j g = a/\a\ respectively) as a function of the proto- 
planet's mass for the four different sizes of the accreted planetesimals of Fig. [6] (corresponding to the dashed line in Fig. ([6])). The 
peaks in the migration timescale correspond to the changes in the sense of migration of the protoplanet. From the two bottom pan- 
els, it can be seen that when the protoplanet's mass is > 10M ffl both timescales become comparable and eventually the migration 
timescale becomes shorter than the accretion timescale, leading to a fast migration of the protoplanet with very little accretion. In 
the case of 1 km planetesimals the protoplanet growth is too slow to gain enough mass to change the type of migration before it is 
lost in the star (when its mass is ~ 30M ffi ). In the case of 0.1 km planetesimals, the protoplanet can grow more massive because its 
accretion rate is larger than in the previous case, then it is able to open a gap in the disc and migration switches to type II, preventing 
it to fall in the star. 



ful. Finally, if they grow up to around Neptune mass while the 
disc is still young but they do not manage to grow fast enough 
their migration accelerates and they end in the central star. 



4.4. Exploring the parameter space of initial conditions 

One may wonder if the examples presented in the previous sec- 
tion are representative of a general case. In order to answer 
this question we apply our planet formation model for a vari- 
ety of protoplanetary discs (including different lifetimes, masses, 
metallicities, etc.) and initial locations for the embryo. 

In this framework, we performed sets of more than 10 000 
simulations. Each set considers a different size for the accreted 
planetesimals (r m = 100, 10, 1 and 0.1 km). For these simu- 



lations we adopt the disc models described in Sec. 2.1.2 The 



initial mass of the embryo is, in all the cases, 0.01 M ffi . The ini- 
tial location of the embryo is varied between 0.2 and 30 AU. 
For an embryo to start at a certain location, we check that the 
initial mass of solids at its location is greater than the mass of 



the embryo. We subtract the mass of the embryo from the initial 
mass of solids of its feeding zone, which means we are assuming 
instantaneous, in situ formation for it. Planets grow by accret- 
ing solids and gas, and can migrate in the disc, according to the 
model we presented in the previous sections. The evolution of 
the eccentricity and inclination of planetesimals are calculated 
solving the differential equations (Eqs. pT| ) - ( [32} ) all through- 
out the disc for every time step. As initial conditions for e and i 
we assume that they are given by Eqs. (58 1 - (59 1, what we have 
called before "an initially cold disc". 



The parameter space explored for the initial mass and life- 
time of protoplanetary discs is schematically shown in Fig. [9] 
The rectangle surrounded by a black box shows the whole set of 
initial conditions studied. Note however that the distribution of 
these parameters is likely not uniform, and not all combinations 
have the same probability of occurrence: long lived discs, as well 
as very low and very massive discs are unlikely. We found that 
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giant planet^] do not form for the initial conditions correspond- 
ing to the grey region in Fig. [9] In fact, for a population of 100 
km and 10 km planetesimals, giant planets do not form under 
any of the initial conditions considered in this work. In the case 
of 100 km population of planetesimals this is due to the fact that 
the accretion time to form a massive solid core is longer than 
the discs' lifetimes. For 10 km planetesimals, some giant planets 
would form if migration were not at work (see Sec.[5]l. However, 
according to disc-planet interaction models, planets do migrate 
and the migration timescale turns out to be much shorter than 
the accretion timescale. Planets are lost in the central star before 
they are able to accrete enough solids to trigger the runaway ac- 
cretion of gas (as seen in the previous section, this would allow 
them to grow faster and switch from type I to type II migration, 
therefore preventing their loss in the central star). 

However, the situation reverts when the accreted planetesi- 
mals are smaller. Fig.[9]shows that massive discs favour the for- 
mation of giant planets. The coloured regions depict the char- 
acteristics of discs where, for a particular initial position, the 
embryo succeeded in growing to a giant planet. This does not 
mean that a giant planet will form in any location of the disc, but 
that formation is possible at certain locations. The yellow region 
shows the parameter space of giant planets that succeeded to sur- 
vive in the disc accreting 1 km planetesimals. Clearly, there is a 
dependence on the discs lifetime: less massive but longer lived 
discs favour the formation of giant planets. In red is plotted the 
situation for the case of r m = 0.1 km accreted planetesimals. As 
smaller planetesimals are accreted faster and more efficiently, the 
disc parameter space that leads to the formation of giant planets 
is bigger. Interestingly, in this case there is no dependence on the 
lifetime of the disc. When the amount of solids present is large 
enough, accretion to form a core with a critical mass proceeds so 
fast that it is always shorter than the lifetime of the discs consid- 
ered here. If the planetesimal disc is massive, accretion can be 
fast enough to enable protoplanets to reach the cross-over mass. 
Massive protoplanets switch from fast type I migration to a much 
slower type II, therefore decelerating and eventually braking be- 
fore reaching the central star, preventing this way the loss of 
planets. 

Very small mass planets are the most abundant outcome in 
all simulations (regardless the size of the accreted planetesimals) 
and their final location could be anywhere in the disc. Planets 
more massive than 10 M®, in general, start their formation be- 
yond the ice line (where solids are more abundant and feeding 
zones are larger), and due to migration they reach the inner re- 
gions of the disc. The final location of these planets extends be- 
tween the inner edge of the disc and 10 AU. Planets in the range 
of 10 2 to 10 3 M ffl are the less abundant ones, and their final lo- 
cations are very different from their initial emplacements. These 
planets undergo a lot of migration, and those that remain are the 
ones that were able to accrete gas fast enough to enter the run- 
away accretion of gas which prevented their loss in the central 
star. Planets in the mass range 10 - 10 2 M ffi are those that undergo 
the largest net displacement from their original location. Most of 
the surviving planets in this mass range have masses lower than 
20 M ffl . These planets, that were not able to reach their cross- 
over masses while they were migrating towards the central star 
are preserved because the disc dissipated before they could fall 
into the star. In the case of r m — 0.1 km, ~ 18% of the simula- 
tions lead to lost planets, 85% of which were not able to reach 



6 By giant planets we mean protoplanets with masses larger than the 
cross-over mass and which survive from disc-planet angular momentum 
exchange and migration. 
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Fig. 9. Lifetime vs. initial solid mass of the protoplanetary discs 
considered in our simulations. Note that the total mass of the 
disc relates to the solid mass of the disc through the gas-dust 
ratio. The rectangle surrounded by a black box schematically 
shows the parameter space explored. The grey rectangle depicts 
the region where giant planets do not form at all. The region of 
the parameter space where protoplanets reached the cross-over 
mass are shown in yellow for protoplanets accreting 1 km plan- 
etesimals, and in red for protoplanets accreting 0. 1 km planetesi- 
mals. For bigger accreted planetesimal (100 and 10 km) no giant 
planets form. 



their cross-over mass. This confirms our analysis of the previ- 
ous section. Planets that are massive enough to undergo rapid 
inward type I migration (> 10 M e ) but whose growth rate is still 
dominated by the accretion of solids are likely to be lost in the 
star. 

In the previous sections we have noticed that results of gi- 
ant planet formation depend upon the dynamical model adopted 
to describe planetesimals' dynamics. The equilibrium approach 
has the advantage of being numerically not very time consum- 
ing. However, it can lead to very different results when compared 
to the case of explicitly solving the differential equations for e 
and i. When solving the equations, we also have the problem 
of setting the initial conditions, which are unknown. To test the 
importance of these assumptions, we performed 10000 simula- 
tions under three different conditions. Each calculation differs 
from the other only in the treatment of the planetesimals e and i. 
This means that for a given set of initial conditions, we calculate 
the formation of the planet three times: one assuming the full 
equilibrium situation for planetesimals, another solving the ex- 
plicit differential equations using as initial conditions for e and 
i the equilibrium value when the stirring timescale of the em- 
bryo equals the nebular gas drag timescale (we call this scenario 
"hot disc") and finally also solving the differential equations but 
using as initial conditions the equilibrium value between mutual 
planetesimal stirring and gas drag damping (we call this scenario 
"cold disc"). 

shows the fraction of surviving planets with respect 
amount of planets formed (surviving planets plus 



10 



tota 



Fig 
to the 

planets lost in the central star). Accreted planetesimals are 0.1 
km size. To make this plot we classified the planets according 
to their final mass in five mass bins. Clearly, in all the cases, 
the most abundant planets are the low mass ones (< 1 M ffi ). The 
equilibrium scenario being the one corresponding to the slowest 
accretion rate, is the one with more planets in the lowest mass 
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5. Discussion 
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Fig. 10. Fraction of surviving planets with final mass in a cer- 
tain mass interval. N tota i is the sum of all the surviving planets 
(nsurv) plus the planets lost in the central star. Final masses of 
the formed planets are separated in five mass bins (M<1 M ffl , 
l<M[M e ]<10, 10<M[M e ]<100, 100<M[M e ]<1000, M>1000 
Me). The solid component of the planets grow by the accretion 
of 0. 1 km planetesimals. Three different approaches in the treat- 
ment of the evolution of the eccentricities and inclinations of 
planetesimals are considered: equilibrium, and explicit calcula- 
tion of the differential equations for a "hot" and "cold" initial 
conditions. 



bin (< 1 M ffi ). On the other extreme, in the case of the cold initial 
disc, accretion of solids at the beginning is more efficient, em- 
bryos grow bigger in a shorter time which, in turn, gives them 
more chances to continue accreting. That is why in the case of 
an initial cold disc planets are more massive. This is evident in 
the histogram: for a cold initial planetesimal disc, the fraction 
of planets in each mass bin is the highest (except, of course, in 
lowest mass bin). The amount of planets in the interval of 10 to 

10 2 M e represents around 2% of the surviving planets. Most of 
them are in the mass range of 10 to 20 M ffl , and just a few in 
the range of 20 to 50 M®. There are no planets in the range of 
50 to 10 2 Me, evidencing the dramatic effect of type I migration 
on these intermediate mass planets. In the mass interval 10 2 - 

10 3 M ffi it can be noticed a profound decrease in the number of 
planets, independently of the approach used for the planetesimal 
dynamics. This can be understood as follows: protoplanets with 
masses greater than 10 2 M ffl are usually in the runaway phase 
of gas accretion. This means that hundreds of Earth masses can 
be accreted in a very short time. Therefore planets in the run- 
away phase easily grow more massive than Jupiter. Only if the 
disc dissipates during this process of accretion, final planetary 
masses can be between Saturn and a few Jupiter masses. To a 
less extent, planet migration has also some consequences in de- 
pleting this bin: although the mass bin in which planet migration 
is most effective in eliminating planets is the 10 - 10 2 Ms range, 
it is still very important in this mass regime, eliminating around 
25% of these planets (especially those closer to the 10 2 M ffl edge 
of the mass interval). 



The calculations presented in this paper focus on the formation 
of planets and how this process is highly affected by the accre- 
tion rate of solids. The accretion rate of solids introduced in- 
tends to be realistic while computationally tractable. Some other 
aspects of the model are therefore rather simplified here; for ex- 
ample we consider the formation of a single planet. The applica- 
tion of these calculations to the formation of planetary systems 
will be presented elsewhere (Alibert et al. 2012). We introduced 
a semi-analytical description for the eccentricities and inclina- 
tions of planetesimals. In this work we calculated explicitly plan- 
etesimals' eccentricities and inclinations, taking into account the 
stirring of the growing planets, the gas drag from the nebula and 
the mutual stirring of the planetesimals themselves. The stirring 
produced by the growing planet excites planetesimals and makes 
their accretion more difficult as it grows. On the other hand, gas 
drag counteracts the stirring effect, an effect more important for 
smaller planetesimals. We have considered three approaches to 
determine the rms eccentricity and inclination of planetesimals: 

- "analytical equilibrium" calculation, as described by Eqs. 

- out of equilibrium, by solving the time evolution of the 
rms eccentricity and inclination of planetesimals (Eqs. (31 

starting from a "cold" planetesimal disc (their ex 



(32 1) 



citation state being the results of planetesimal-planetesimal 
interaction and gas drag only), and 
- out of equilibrium, starting from a "hot" disc, where plan- 
etesimals are already excited by the planetary embryo. 

We have shown that these three approaches lead to different 
accretion rates, the difference depending on the planetesimal 
size, and being more important as a result of the migration and 
disc evolution feedbacks. For large planetesimals, the three ap- 
proaches lead to similar results, but the accretion rate of solids 
under this hypothesis is very small, preventing the formation 
of massive planets. On the other hand, for low mass planetes- 
imals, the excitation state of planetesimals as derived from the 
"analytical equilibrium" calculations, and the excitation state of 
planetesimals following the second or third approach, even after 
a long time (when the equilibrium solution of Eqs. ( [31} - ( [32} 
is reached), are different, resulting in different accretion rate of 
solids. In particular, the ratio of the eccentricity to the inclination 
can deviate substantially from the 1/2 ratio assumed in the first 
approach. 

As a consequence of the size dependence of gas drag, small 
planetesimals are easier to accrete, leading to a faster formation 
of planets. We have shown that in the framework of the model 
hypothesis outlined in Sect. [2] the formation of planets rang- 
ing from a fraction of an Earth mass to several Jupiter masses 
can only be accomplished under the assumption of a population 
of small planetesimals or, more generally, if planetesimals are, 
by any process, maintained in a "cold" dynamical state. Similar 
conclusions were already identified in other works (e.g. Fortier et 
al. 2007, 2009, Benvenuto et al. 2009), where it was shown that, 
in order to achieve the formation of the giant planets of the so- 
lar system in less than 10 Myr most of the accreted planetesimals 
have to be small. Those models, however, assumed in situ forma- 
tion, and did not consider a consistent calculation of the planet's 
final mass (the computations were stopped when the masses of 
the giant planets of the solar system were reached). Although the 
approach for the calculation of planet formation (regarding the 
accretion of gas and solids) is similar, our work accounts for the 
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Fig. 11. In situ formation. Fraction of planets with final mass in a certain mass interval. N tota i represents the total amount of 
planets formed, independently of their mass. Planets are separated in five mass bins (M<1 M e , l<M[M e ]<10, 10<M[M e ]<100, 
100<M[M e ]<1000, M>1000 M ffi ). Accreted planetesimals have a radius of 100 km (red), 10 km (green), 1 km (blue) and 0.1 km 
(magenta). Most of the final masses are in the first mass bin. The inset caption represents the magnified histogram for the other mass 
bins. 



migration of the planets and their final masses are determined by 
the coupled evolution of the disc and the planet. 

In order to compare our results with the afore-mentioned 
works, we performed simulations where migration is switched 
off. We calculate the in situ formation of 10000 planets that ac- 
crete planetesimals of 100, 10, 1 and 0.1 km in radius. Indeed, 
as can be seen from Fig.[TT[ small planets are the most abundant 
ones. Larger planets form when there are small planetesimals 
to accrete (smaller than r m = 10 km). Planets in the mass range 
10 - 10 2 M ffl (which growth is not yet dominated by the accretion 
of gas) are produced in about the same fraction as giant planets 
in the mass range > 10 3 M®, which growth is dominated by the a 
runaway accretion of gas. The runaway gas accretion leads to the 
difficulty of forming planets in the mass range of 10 2 - 10 3 M ffi : 
only a fine tuned timing effect, the gas disc disappearing when 
the planet is in this mass range, can lead to planets in the Saturn- 
Jupiter domain. 

If we focus on the mass bin 10 to 10 2 M e , one difference 
between the migration and in situ scenario is that, while in the 
former there are no planets in the mass range between 50 and 
100 M ffi (because they are all engulfed by the central star), in the 
in situ case planets are distributed in the whole mass spectrum 
of the bin (however, most of them are in the lowest mass range 
[10 to 20 M ffl ]). It is also interesting to compare for these two 



scenarios the mass percentage of solids for planets between 5 - 
10 2 IVLj. As it can be seen in table [3] planets that migrate have a 
more massive core than planets formed in situ. In general, when 
planets migrate they have access to regions of the disc that have 
not been depleted of planetesimals. Therefore, the solids surface 
density is higher and so is the accretion rate, which results in the 
formation of a more massive core. 

In this paper we have considered that the opacity of the 
envelope corresponds to the full interstellar medium opacity. 
However, several works (e.g. Pollack et al. 1996, Podolak 2003, 
Hubickyj et al. 2005, Movshovitz & Podolak 2008) suggest that 
opacity in the planet's envelope should be much smaller, lead- 
ing to a faster formation. As a complementary result, we com- 
puted one in situ population and one population including mi- 
gration (10 000 planets each) reducing the opacity to 2% of that 
of the interstellar medium. We found that with a reduced opac- 
ity the loss of planets decreases, increasing in every mass bin 
the number of planets that survive in the protoplanetary disc. 
Nevertheless, the shape of the mass distribution is very simi- 
lar to that shown in Fig. 10 the fraction of surviving planets 



decreases with the mass of the planet, with a minimum in the 
interval 10 2 - 10 3 M e (although with the reduced opacity the 
fraction of planets in this bin is four times larger than with the 
full opacity), and raising again in the last mass interval. In ta- 
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Table 3. For planets whose total mass is between 5 and 10 2 M®, 
this table shows the percentage of the total mass contained in the 
solid core. A comparison between a migration and non migration 
scenario is also made. In the migration case there are no planets 
with masses greater than 50 M©. 



Mass Interval 


M corc /M tota i 


x 100 


[M e ] 


With Migration 


In situ 


5-10 


98 


94 


10-20 


95 


86 


20-30 


88 


75 


30-40 


78 


63 


40-50 


72 


56 


50-60 




56 


60-70 




43 


70-80 




37 


80-90 




34 


90 - 100 




39 



Table 4. The same as in |3jbut considering that the planet's en- 
velope opacity is only 2% of the interstellar medium opacity. 



Mass Interval 


M corc /M total 


x 100 


[M e ] 


With Migration 


In situ 


5-10 


96 


80 


10-20 


90 


72 


20-30 


79 


62 


30-40 


71 


58 


40-50 


68 


46 


50-60 


50 


44 


60-70 


51 


34 


70-80 


48 


29 


80-90 




32 


90 - 100 


32 


20 



ble|4]we show the solids mass fraction for planets with masses 
between 5 to 100 M®, considering both in situ and migration sce- 
nario. In both cases, the fraction of solids is smaller than in the 
corresponding case calculated with full opacity. 

Using the full model (including migration and disc evolu- 
tion), giant planet formation by accretion of 100 km planetes- 
imals is quite unlikely, if not impossible. In our simulations, 
to actually form giant planets we had to reduce the planetesi- 
mal size to the 0.1 km - 1 km radius range. Such a conclusion 
raises the question of the most likely planetesimal size during the 
epoch of planet formation. Recent studies on planetesimal for- 
mation however give different conclusions. On one hand, models 
that explain the formation of planetesimals by direct collapse in 
vortices in turbulent regions (e.g. Johansen et al. 2007) predict 
a fast formation of very big planetesimals (r,„ > 100 km). We 
note however that this formation process may not be totally effi- 
cient and only a small fraction of solids initially present in pro- 
toplanetary discs is likely to end up in such big planetesimals. 
The conclusions of Johansen et al. (2007) are consistent with the 
results of Morbidelli et al. (2009) on the initial function of plan- 
etesimals. On the other hand, a recent study of Windmark et al. 
(2012) shows that direct growth of planetesimals via dust colli- 
sions can lead to the growth of 0.1 km planetesimals. In addi- 
tion, initially small planetesimals show better matches to the ob- 
served size distribution of objects in the asteroid belt and among 
the TNOs: Weidenschilling (2011) shows that the size distribu- 
tion currently observed in the asteroid belt in the range of 10 
to 100 km can be better explained by an initial population of 
0.1 km planetesimals. Kenyon and Bromley (2012) conclude, by 
combining observations of the hot and cold populations of TNOs 



with time constraints on their formation process, that TNOs form 
from a massive disc mainly composed of 1 km planetesimals. 
More investigations on the formation of planetesimals, and plan- 
etary embryos, are definitely required in order to test the viability 
of planetary core formation by accretion of low mass planetesi- 
mals. We note finally that what is important in the work we have 
presented now is not the initial mass function of planetesimals, 
but their mass function at the time of planet formation. The two 
quantities are likely to differ due to planetesimal-planetesimal 
collisions, and the resulting mass growth and/or fragmentation. 

6. Conclusions 

In this paper we presented calculations of planetary formation 
considering the formation of a single planet at each time, and 
starting with embryos of 0.01 M ffl . Our simulations consist in the 
calculation of the formation of a planet, including its growth in 
mass by accretion of solids and gas, its migration in the disc, 
and the evolution of the disc until the gas component of the 
disc is dispersed. When the nebular gas is gone, simulations are 
stopped. Therefore, the subsequent growth of the planets, by ac- 
cretion of residual planetesimals or collisions among embryos 
if we were considering planetary systems, is not considered. 
During their formation, the growth of the planets is calculated 
self-consistently coupling the accretion of solids and gas. The 
accretion of solids is computed assuming the particle-in-a-box 
approximation, and computing the excitation state of planetesi- 
mals, which in turn regulates to a great extent the accretion rate 
of the planet. The accretion of gas is computed solving the dif- 
ferential equations that govern the evolution of the structure of 
the planet. Finally, protoplanets grow in an evolving protoplan- 
etary disc, which density, temperature and pressure is calculated 
at every time step. 

The combination of oligarchic growth (for the solid compo- 
nent of the planet) with the migration of the planet has severe 
consequences for protoplanets that are able to grow up to a few 
tens Earth masses: these planets tend to collide with the central 
star (or at least to migrate to the innermost location of the pro- 
toplanetary disc). Indeed, planets that are between 10 and 100 
Earth masses are usually undergoing very rapid inward type I 
migration, but are not massive enough to switch to a slower, type 
II migration. In our simulations, the only surviving planets in the 
range of 10 to 20 Earth masses correspond to cases where the gas 
component of the disc dissipates during their growth, preventing 
them to fall into the star. On the opposite, if the solid core grows 
fast enough, it enables the accretion of large amounts of gas, 
when the critical mass is reached. At this point, the runaway of 
gas ensures an extremely quick growth in the mass of the planet, 
and the planet migration rate decreases. 

In the model we have presented in this paper, we have as- 
sumed a uniform population of small planetesimals which size 
remains unchanged during the whole formation of the planet. 
Most probably the initial population of planetesimals in proto- 
planetary discs is not uniform in size, but follows a size distri- 
bution. We have shown, however, that without small planetesi- 
mals giant planet formation is difficult to explain, at least in the 
way we understand it now. However, even with an initial pop- 
ulation of small planetesimals, the collisions among themselves 
are likely to be disruptive as soon as their random velocities start 
to be excited by a planetary embryo. Therefore, it is also unlikely 
that an initial population of only small planetesimals can be used 
to explain the formation of giant planets. Moreover, even assum- 
ing this to be true, only a few amount of planets in the range of 
several tens Earth masses to a few Jupiter masses can be formed. 
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In addition, small mass planetesimals are subject to large ra- 
dial drift, as a result of gas drag. Planetesimal drift can have 
positive or negative consequences in the formation of plane- 
tary systems, as has been shown by Guilera et al. (2010, 2011). 
Similarly, fragmentation and coagulation can hasten or delay 
planet formation as a whole: fragments of smaller mass are eas- 
ier to accrete but, if they also can leave a planet's feeding zone 
very quickly as a result of gas drag. Finally, in a planetary sys- 
tem, fragments that are not accreted by the embryo that gen- 
erated them can be accreted by another embryo located in an 
innermost region. It is not clear what the possible outcomes of 
putting together planetesimals drifting, fragmentation and many 
embryos forming in an evolving disc could be. This however 
represents a very important step in the understanding of the first 
stage of planet formation. 

Because most of the accretion of solids should, at some 
point, be dominated by small planetesimals or fragments, our 
calculations can be understood as a description of that stage. 
An interesting scenario to analyse, in particular if planetesi- 
mals are born massive (Johansen et al. 2007, Morbidelli et al. 
2009) would be the following. An initial population dominated 
by ~ 100 km planetesimals would prevent a fast growth of the 
embryos at the beginning, a time during which the planetary 
embryos would suffer only of a little migration and the proto- 
planetary disc would evolve, reducing its gas surface density. 
Fragments (smaller planetesimals), resulting from collisions be- 
tween big planetesimals, would start to be generated later (the 
timescale of fragmentation of 100 km planetesimals affected by 
the stirring of Moon to Mars mass embryos is of the order 10 6 
years), therefore accelerating the formation of the embryo in a 
later stage of the disc evolution. The collisional cascade would 
probably still produce small fragments fast enough to help the 
growth of an embryo, even if they leave the feeding zone very 
fast. Therefore, protoplanets could grow by the accretion of frag- 
ments, not necessarily generated by themselves, but generated 
by another distant embryo. Putting together these different pro- 
cesses will give us a better insight on the formation process and 
would help us to constrain, from planetary formation models, 
possible initial size distributions for planetesimals. 

It is also important to mention that we have not considered 
the possibility of planetesimal-driven migration. Although in 
general planetesimal-driven migration acts on a longer timescale 
that Type I migration, recently Ormel et al. (2012) find that 
planetesimal-driven migration can have a mild effect on mid- 
sized planets in massive planetesimal discs, competing with 
Type I migration. 

The main conclusions of our work are that formation of giant 
planets in the framework of the sequential accretion model needs 
the presence of unexcited planetesimals. One obvious way to de- 
excite planetesimals is through gas drag, but this requires this 
latter to be efficient, which in turn translates to low mass plan- 
etesimals. These planetesimals can be primordial or fragments 
of originally bigger ones. But, at some point, small boulders are 
needed to build protoplanetary massive cores before the dissi- 
pation of the disc. The combination of migration and oligarchic 
growth, on the other hand, prevents the formation of intermedi- 
ate mass planets. However, this result can change when consid- 
ering the formation of planets in planetary systems where their 
gravitational interactions are taken into account. Captures in res- 
onances can prevent planets from colliding with the central star, 
preserving them in planetary systems. Exploring this effects will 
be the subject of future works. 
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